MGP Calculation (rCTP derivation)

FALSE -- Attaching packages --------------------------------------- tidyverse 1.3.0 --

FALSE v ggplot2 3.3.2     v purrr   0.3.4
FALSE v tibble  3.0.3     v dplyr   1.0.2
FALSE v tidyr   1.1.2     v stringr 1.4.0
FALSE v readr   1.4.0     v forcats 0.5.0

FALSE -- Conflicts ------------------------------------------ tidyverse_conflicts() --
FALSE x tidyr::extract()   masks magrittr::extract()
FALSE x dplyr::filter()    masks stats::filter()
FALSE x dplyr::lag()       masks stats::lag()
FALSE x purrr::set_names() masks magrittr::set_names()

FALSE Loading required package: Hmisc

FALSE Loading required package: lattice

FALSE Loading required package: survival

FALSE Loading required package: Formula

FALSE 
FALSE Attaching package: 'Hmisc'

FALSE The following objects are masked from 'package:dplyr':
FALSE 
FALSE     src, summarize

FALSE The following objects are masked from 'package:base':
FALSE 
FALSE     format.pval, units

FALSE Loading required package: SparseM

FALSE 
FALSE Attaching package: 'SparseM'

FALSE The following object is masked from 'package:base':
FALSE 
FALSE     backsolve

Now that we’ve loaded in the necessary libraries we can set up the marker gene lists to run the MGP algorithm. These marker gene lists will be used to describe each of the cell types in the relative cell-type proportions (rCTPs) we are calculating for each cohort.

There are marker lists calculated from the cingulate gyrus (CgG) at the resolution of inhibitory/excitatory (or GABAergic/glutamatergic) neuronal cells, marker lists calculated from the CgG at the resolution of neuronal subclasses (SST cells, IT cells) and marker lists calculated from middle temporal gyrus (MTG) at the resolution of neuronal subclasses. These lists are denoted as “in_ex_CgG”, “subclass_CgG”, “subclass_MTG” respectively. Additionally there are validated published marker gene lists for comparison from the Lake et. al (2017) and Darmanis et al. (2015) datasets (defined at the inhibitory/excitatory resolution), hereafter referred to as “lake” and “darmanis”.

marker_lists <- c("in_ex_CgG","subclass_CgG", "subclass_MTG", "lake", "darmanis")
for(marker in marker_lists){
  setwd('./markerLists')
  markers_df <-  readRDS(paste0(marker, '.rds'))
  assign(marker, markers_df)
  setwd('../')
}

Now that we’ve loaded in all the marker lists we want to load in all our cohorts count matrices. The cohorts were defined in pre-processing and we will continue to use those definitions: the Religious Orders Study and Rush Memory and Aging Project cohort (delineated ROSMAP from here on) which has dorsolateral prefrontal cortex samples, the Mayo Clinic cohort, hereafter referred to as MAYO, which consists of temporal cortex samples and the Mount Sinai Brain Bank cohort (referred to as MSBB), which has expression data sampled from Brodmann area 10, 22, 36 and 44 (written BM10, BM22, BM36, BM44 hereafter).

#load in new cohorts QC-ed data, final count matrices and convert them to dataframes with HGNC symbols instead of ENSEMBL ids
cohorts <- c("ROSMAP", "MAYO", "MSBBM10", "MSBBM22", "MSBBM36", "MSBBM44")
for(cohort in cohorts){
  if(str_detect(cohort, "MSBB")){
    cohort <- gsub('MSB', '', cohort) 
  }
  print(cohort)
  matrix_name <- paste0(cohort, "_matrix")
  
  filename <- paste0(matrix_name, ".rds")
  setwd('./finalCountMatrices')
  matrix <- readRDS(filename)
  setwd('../')
  assign(matrix_name, matrix)

  count_df <- as.data.frame(matrix)
  count_df <- tibble:: rownames_to_column(count_df, var="Gene")
  setwd('./geneAnno')
  gene_anno <- readRDS("gene_anno.rds")
  setwd('../')
  final_df <- merge(count_df, gene_anno)
  new_gene <- final_df$Hgnc_Gene %>% make.names(unique = T) #this is ONE way of dealing with the duplicates, just making them into separate, unique names
  final_df$new_gene <- new_gene
  n_col <- ncol(count_df)
  mgp_df <- final_df[,c(n_col+2,2:n_col)]
  count_df <- mgp_df %>% rename(Gene = new_gene)
  df_name <- paste0(cohort, "_count_df")
  assign(df_name, count_df)
}

## [1] "ROSMAP"
## [1] "MAYO"
## [1] "BM10"
## [1] "BM22"
## [1] "BM36"
## [1] "BM44"

Now that we have our marker lists and our cohort data loaded and converted to HGNC symbols, we want to get the list of common marker genes that are found within each cohort. Each marker gene list defines a cell type through a list of marker genes that are highly expressed and uniquely expressed in that cell type (to an extent). However, not all of the marker genes that define cell type proportions are found in each of the cohorts, for example: geneX is used to define celltypeX, and is found in cohorts BM10, BM22, and ROSMAP. Since geneX is not found in MAYO, BM36 or BM44 we do not want to include it in the list of marker genes used to define celltypeX as these relative cell type proportions (rCTPs) will then be less comparable across cohorts. As such we’re going to run a QC function for the markerGeneProfile (MGP) method on all of the markers and then take all the markers found across all cohorts, and drop all the celltypes that have 4 or less genes defining them.

Note that we will add a “new” cell type marker list “subclass_both” just to run the MGP algorithm. This is not really a new marker list, it simply applies the subclass_MTG and subclass_CgG lists to the 6 cohorts based on their brain region specificity. This means the marker list “subclass_both” will result in the MGP algorithm running with the subclass_MTG marker list on the temporal cortex samples (MAYO, BM36, BM22) because the MTG marker list is derived from brain regions nearest to the temporal cortex, and running the subclass_CgG marker list on the prefrontal cortex samples (ROSMAP, BM10, BM44) as that marker list is derived from the region closest to prefrontal. We will then take the common markers found in the temporal cortex samples (MAYO, BM36, BM22) in the MTG list to run the MGP algorithm and then take the common markers found in the prefrontal cortex samples (ROSMAP, BM10, BM44) in the CgG list to run the MGP algorithm.

#we define our QC algorithm such that it returns a dataframe with the cell type, markers_used (list of marker genes used ' per cell type), removed_marker_ratios (list of removed marker ratios per cell type) and percent_variance_PC1 (list of variance explained by the first PC per cell type)
mgpQCMetrics <-function(count_df, mgp_markers){
  
  if(!all(c("Gene") %in% colnames(count_df))){
    stop("The count_df argument must be a df with a column named Gene (HGNC gene symbols)")
  }
  mgp_est<- markerGeneProfile::mgpEstimate(exprData=count_df,
                                                  genes=mgp_markers,
                                                  geneColName="Gene",
                                                  outlierSampleRemove=F, # should outlier samples removed. This is done using boxplot stats.
                                                  geneTransform =NULL, #function(x){homologene::mouse2human(x)$humanGene}, # this is the default option for geneTransform
                                                  groups=NULL, #if there are experimental groups provide them here. if not desired set to NULL
                                                  seekConsensus = FALSE, # ensures gene rotations are positive in both of the groups
                                                  removeMinority = TRUE)
  i= 0
  for(cell in names(mgp_markers)){
    i = i + 1
    cells_df <- mgp_est$usedMarkerExpression[i] %>% as.data.frame()
    masterlist <- paste0(rownames(cells_df), collapse=', ')
    num_markers <- length(rownames(cells_df))
    rm_marker_ratios <- mgp_est$removedMarkerRatios[i]
    if(!is.null(mgp_est$trimmedPCAs[[i]])){
      percent_variance <- ((summary(mgp_est$trimmedPCAs[[i]]))[6]) %>% as.data.frame()
      percent_variance_PC1 <- percent_variance[2,1]
    }
    else{
      percent_variance_PC1 <- NA
    }
    if(i==1){
      master_df <- data.frame( "markers_used" = masterlist, 
                        "removed_marker_ratios" = rm_marker_ratios,
                        "percent_variance_PC1" = percent_variance_PC1, 
                        "num_markers" = num_markers)  
    }
    else{
      df <- data.frame( "markers_used" = masterlist, 
                        "removed_marker_ratios" = rm_marker_ratios,
                        "percent_variance_PC1" = percent_variance_PC1, 
                        "num_markers" = num_markers)
      master_df <- rbind(master_df, df)
    }
  }
  master_df <- tibble::rownames_to_column(master_df, var = "celltype")
  return(master_df)
}


#calculate QC metrics for each marker gene list for each cohort
cohorts <- c("ROSMAP", "MAYO", "MSBBBM10", "MSBBBM22", "MSBBBM36", "MSBBBM44")
marker_lists <- c("in_ex_CgG","subclass_CgG", "subclass_MTG", "subclass_both", "lake", "darmanis")

for(markers in marker_lists){
  for(cohort in cohorts){
    if(str_detect(cohort, "MSBB")){
      cohort <- gsub('MSBB', '', cohort) 
    }
    print(cohort)
    df_name <- paste0(cohort, "_count_df")
    if(markers != "subclass_both"){
       mgpResult <- mgpQCMetrics(get(df_name), mgp_markers = get(markers))
    }
    else{
      #if subclass_both specified run CgG markers on prefrontal and MTG on temporal
      if(cohort == "ROSMAP" | cohort =="BM10" | cohort =="BM44"){
        mgpResult <- mgpQCMetrics(get(df_name), mgp_markers = subclass_CgG)
      }
      else{
         mgpResult <- mgpQCMetrics(get(df_name), mgp_markers = subclass_MTG)
      }
    }
    mgpResult$cohort <- cohort
    mgpName <- paste0("mgpQCResults",cohort, markers)
    assign(mgpName, mgpResult)
    setwd(paste0('./mgpQC_',markers))
    saveRDS(get(mgpName), paste0(mgpName, ".rds"))
    setwd('../')
    print(mgpName)
    assign(mgpName, mgpResult)
    if(cohort =="ROSMAP"){
      all_cohorts_QC <- mgpResult
    }
    else{
      all_cohorts_QC <-rbind(all_cohorts_QC, mgpResult)
    }
  }
  all_cohorts_QC$cohort <- factor(all_cohorts_QC$cohort, levels = c("ROSMAP", "BM10", "BM44", "BM22", "BM36", "MAYO"))
  all_cohorts_QC <- arrange(all_cohorts_QC, cohort)
  setwd(paste0('./mgpQC_',markers))
  saveRDS(all_cohorts_QC, "all_cohorts_QC.rds")
  setwd('../')
}

## [1] "ROSMAP"
## [1] "mgpQCResultsROSMAPin_ex_CgG"
## [1] "MAYO"
## [1] "mgpQCResultsMAYOin_ex_CgG"
## [1] "BM10"
## [1] "mgpQCResultsBM10in_ex_CgG"
## [1] "BM22"
## [1] "mgpQCResultsBM22in_ex_CgG"
## [1] "BM36"
## [1] "mgpQCResultsBM36in_ex_CgG"
## [1] "BM44"
## [1] "mgpQCResultsBM44in_ex_CgG"
## [1] "ROSMAP"
## [1] "mgpQCResultsROSMAPsubclass_CgG"
## [1] "MAYO"
## [1] "mgpQCResultsMAYOsubclass_CgG"
## [1] "BM10"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: L6b

## [1] "mgpQCResultsBM10subclass_CgG"
## [1] "BM22"
## [1] "mgpQCResultsBM22subclass_CgG"
## [1] "BM36"
## [1] "mgpQCResultsBM36subclass_CgG"
## [1] "BM44"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: L6b

## [1] "mgpQCResultsBM44subclass_CgG"
## [1] "ROSMAP"
## [1] "mgpQCResultsROSMAPsubclass_MTG"
## [1] "MAYO"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: L5/6 NP, L4 IT

## [1] "mgpQCResultsMAYOsubclass_MTG"
## [1] "BM10"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: L5/6 NP

## [1] "mgpQCResultsBM10subclass_MTG"
## [1] "BM22"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: Microglia, L4 IT

## [1] "mgpQCResultsBM22subclass_MTG"
## [1] "BM36"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: L5/6 NP, SST

## [1] "mgpQCResultsBM36subclass_MTG"
## [1] "BM44"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: L5/6 NP

## [1] "mgpQCResultsBM44subclass_MTG"
## [1] "ROSMAP"
## [1] "mgpQCResultsROSMAPsubclass_both"
## [1] "MAYO"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: L5/6 NP, L4 IT

## [1] "mgpQCResultsMAYOsubclass_both"
## [1] "BM10"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: L6b

## [1] "mgpQCResultsBM10subclass_both"
## [1] "BM22"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: Microglia, L4 IT

## [1] "mgpQCResultsBM22subclass_both"
## [1] "BM36"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: L5/6 NP, SST

## [1] "mgpQCResultsBM36subclass_both"
## [1] "BM44"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: L6b

## [1] "mgpQCResultsBM44subclass_both"
## [1] "ROSMAP"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: Ex5, Ex6, Ex8, In2

## [1] "mgpQCResultsROSMAPlake"
## [1] "MAYO"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: Ex4, Ex5, Ex6, Ex8, In4, In6

## [1] "mgpQCResultsMAYOlake"
## [1] "BM10"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: Ex4, In1

## [1] "mgpQCResultsBM10lake"
## [1] "BM22"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: Ex4, Ex5

## [1] "mgpQCResultsBM22lake"
## [1] "BM36"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: Ex2, Ex8, In2

## [1] "mgpQCResultsBM36lake"
## [1] "BM44"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: Ex2, Ex8, In2

## [1] "mgpQCResultsBM44lake"
## [1] "ROSMAP"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: 6

## [1] "mgpQCResultsROSMAPdarmanis"
## [1] "MAYO"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: 6

## [1] "mgpQCResultsMAYOdarmanis"
## [1] "BM10"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: 6

## [1] "mgpQCResultsBM10darmanis"
## [1] "BM22"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: 6

## [1] "mgpQCResultsBM22darmanis"
## [1] "BM36"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: 6

## [1] "mgpQCResultsBM36darmanis"
## [1] "BM44"

## Warning in markerGeneProfile::mgpEstimate(exprData = count_df, genes = mgp_markers, : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: 6

## [1] "mgpQCResultsBM44darmanis"

We now have a wide variety of QC results for each of the marker gene lists and each of the cohorts. Let’s plot some of the data to see what our results look like, before we only take the genes found commonly across each cohort.

FALSE Saving 7 x 5 in image

FALSE Saving 7 x 5 in image

FALSE Saving 7 x 5 in image

FALSE Saving 7 x 5 in image

FALSE Saving 7 x 5 in image

FALSE Saving 7 x 5 in image

These plots clearly show the subclass_both marker list runs the subclass_MTG list on temporal samples and the subclass_CgG list on prefrontal samples.

We can also see we don’t have many markers in the darmanis and lake sets, but let’s get all the common marker genes anyways and plot again — keeping in mind that we need to take the common list within BM22,BM36, and MAYO and ROSMAP, BM10 and BM44 separately for the subclass_both marker list.

FALSE [1] "in_ex_CgG"
FALSE [1] "subclass_CgG"
FALSE [1] "subclass_MTG"
FALSE [1] "subclass_both"
FALSE [1] "lake"
FALSE [1] "darmanis"

Now that we’ve got the common list among the relevant cohorts, let’s plot.

FALSE [1] "in_ex_CgG"

FALSE Saving 7 x 5 in image

FALSE [1] "subclass_CgG"

FALSE Saving 7 x 5 in image

FALSE [1] "subclass_MTG"

FALSE Saving 7 x 5 in image

FALSE [1] "subclass_both_CgG"

FALSE Saving 7 x 5 in image

FALSE [1] "subclass_both_MTG"

FALSE Saving 7 x 5 in image

FALSE [1] "lake"

FALSE Saving 7 x 5 in image

FALSE [1] "darmanis"

FALSE Saving 7 x 5 in image

We lost too many markers in darmanis & lake for those common lists to be useful marker lists, so we’ll stick to using the darmanis and lake marker lists but now use “in_ex_CgG_common_final”, “subclass_CgG_common_final”, “subclass_MTG_common_final”, “subclass_both_MTG_common_final” and “subclass_both_CgG_common_final” as our new marker lists to run the MGP algorithm. This will encourage comparability across cohorts of our results and allow us to perform mega-analysis of the rCTPs across the cohorts/brain regions.

Let’s run the MGP algorithm with our marker lists.

mgpCalc<-function(count_df, meta_df, markers){
  # calculate MGPs per sample
  estimations_human_markers<- mgpEstimate(exprData=count_df,
                                          genes=markers,
                                          geneColName="Gene",
                                          outlierSampleRemove=F, # should outlier samples removed. This is done using boxplot stats.
                                          geneTransform = NULL,
                                          #function(x){homologene::mouse2human(x)$humanGene}, # this is the default option for geneTransform
                                          groups=NULL, #if there are experimental groups provide them here. if not desired set to NULL
                                          seekConsensus = FALSE, # ensures gene rotations are positive in both of the groups
                                          removeMinority = TRUE)
  mgp_info <- list("mgp_df"= count_df, "estimations_human_markers" = estimations_human_markers)
  estimations_human_markers <- mgp_info$estimations_human_markers
  
  # matrix of mgp estimates per cell type, column name SST stores SST MGPs
  mgp_est <- estimations_human_markers$estimates %>% as.data.frame() 
  colnames(mgp_est) <- names(markers)
  mgp_est <- mgp_est %>% tibble::rownames_to_column(var = 'projid')
  
  # merge mgp data frame with sample metadata data frame
  mgp_est = merge(meta_df, mgp_est, by = 'projid')
  return(mgp_est)
}


#cohort can be ROSMAP, MAYO, MSBBBM10, MSBBBM22, MSBBBM26, MSBBBM44
mgpsForAMPAD<-function(cohort, markers){
  if(cohort == "ROSMAP"){
    covars <- c("sex","age_death")
    }
  else{
    covars <- c("msex","AgeAtDeath")
  }
  setwd('./QCpipelineResults')
  v_name <-  paste0("v_", cohort)
  voom <- readRDS(paste0(v_name, ".rds"))
  assign(v_name, voom)
  setwd('../')
  pheno_df <- voom$targets
  if(cohort == "ROSMAP" || cohort == "MAYO"){
    pheno_df<- pheno_df %>% 
      rename(
        projid = SampleID,
      )
  }
  else{
    pheno_df<- pheno_df %>% 
      rename(
        projid = sampleIdentifier,
      )
  }
  count_df_name <- paste0(cohort, "_count_df")
  mgp_result <- mgpCalc(get(count_df_name), pheno_df, markers)
  model.data <- mgp_result
  
  return(list("model" = model.data, "covars" = covars))
}

cohorts <- c("ROSMAP", "MAYO", "BM10", "BM22", "BM36", "BM44")
marker_lists <- c("in_ex_CgG_common_final","subclass_CgG_common_final", "subclass_MTG_common_final","subclass_both", "lake", "darmanis")
for(markers in marker_lists){
  print(markers)
  for(cohort in cohorts){
    if(markers != "subclass_both"){
      mgp_result <- mgpsForAMPAD(cohort, get(markers))
    }
    else{
      #if subclass_both specified run CgG markers on prefrontal and MTG on temporal
      if(cohort == "ROSMAP" | cohort =="BM10" | cohort =="BM44"){
        mgp_result <- mgpsForAMPAD(cohort, subclass_both_CgG_common_final)
      }
      else{
        mgp_result <- mgp_result <- mgpsForAMPAD(cohort, subclass_both_MTG_common_final)
      }
    }
    mgp_name <- paste0("mgp_",cohort)
    print(mgp_name)
    assign(mgp_name, mgp_result)
    folder_name <- gsub('_common_final', '', markers)
    setwd(paste0('./mgpResults_',folder_name))
    saveRDS(get(mgp_name), paste0(mgp_name, ".rds"))
    setwd('../')
  }
}

## [1] "in_ex_CgG_common_final"
## [1] "mgp_ROSMAP"
## [1] "mgp_MAYO"
## [1] "mgp_BM10"
## [1] "mgp_BM22"
## [1] "mgp_BM36"
## [1] "mgp_BM44"
## [1] "subclass_CgG_common_final"
## [1] "mgp_ROSMAP"
## [1] "mgp_MAYO"
## [1] "mgp_BM10"
## [1] "mgp_BM22"
## [1] "mgp_BM36"
## [1] "mgp_BM44"
## [1] "subclass_MTG_common_final"
## [1] "mgp_ROSMAP"
## [1] "mgp_MAYO"
## [1] "mgp_BM10"
## [1] "mgp_BM22"
## [1] "mgp_BM36"
## [1] "mgp_BM44"
## [1] "subclass_both"
## [1] "mgp_ROSMAP"
## [1] "mgp_MAYO"
## [1] "mgp_BM10"
## [1] "mgp_BM22"
## [1] "mgp_BM36"
## [1] "mgp_BM44"
## [1] "lake"

## Warning in mgpEstimate(exprData = count_df, genes = markers, geneColName = "Gene", : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: Ex5, Ex6, Ex8, In2

## [1] "mgp_ROSMAP"

## Warning in mgpEstimate(exprData = count_df, genes = markers, geneColName = "Gene", : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: Ex4, Ex5, Ex6, Ex8, In4, In6

## [1] "mgp_MAYO"

## Warning in mgpEstimate(exprData = count_df, genes = markers, geneColName = "Gene", : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: Ex4, In1

## [1] "mgp_BM10"

## Warning in mgpEstimate(exprData = count_df, genes = markers, geneColName = "Gene", : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: Ex4, Ex5

## [1] "mgp_BM22"

## Warning in mgpEstimate(exprData = count_df, genes = markers, geneColName = "Gene", : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: Ex2, Ex8, In2

## [1] "mgp_BM36"

## Warning in mgpEstimate(exprData = count_df, genes = markers, geneColName = "Gene", : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: Ex2, Ex8, In2

## [1] "mgp_BM44"
## [1] "darmanis"

## Warning in mgpEstimate(exprData = count_df, genes = markers, geneColName = "Gene", : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: 6

## [1] "mgp_ROSMAP"

## Warning in mgpEstimate(exprData = count_df, genes = markers, geneColName = "Gene", : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: 6

## [1] "mgp_MAYO"

## Warning in mgpEstimate(exprData = count_df, genes = markers, geneColName = "Gene", : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: 6

## [1] "mgp_BM10"

## Warning in mgpEstimate(exprData = count_df, genes = markers, geneColName = "Gene", : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: 6

## [1] "mgp_BM22"

## Warning in mgpEstimate(exprData = count_df, genes = markers, geneColName = "Gene", : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: 6

## [1] "mgp_BM36"

## Warning in mgpEstimate(exprData = count_df, genes = markers, geneColName = "Gene", : Several cell types have a high proportion (>0.4) of their genes filtered out. Excercise caution.
## Problematic cell types are: 6

## [1] "mgp_BM44"

Now we’ve calculated MGPs for each of the cohorts and each of the common marker gene lists. We want to Z-score these MGPs so we can compare them in a mega-analysis within each marker gene list.

cohorts <- c("ROSMAP", "MAYO", "BM10", "BM22", "BM36", "BM44")
marker_lists <- c("in_ex_CgG_common_final","subclass_CgG_common_final", "subclass_MTG_common_final","subclass_both", "lake", "darmanis")
#Z score mgps per subject
for(markers in marker_lists){
  folder_name <- gsub('_common_final', '', markers)
  for(cohort in cohorts){
    mgp_name <- paste0("mgp_",cohort)
    setwd(paste0('./mgpResults_',folder_name))
    mgp <- readRDS(paste0(mgp_name, ".rds"))
    mgp_df <- mgp$model
    if(markers != "subclass_both"){
      markers_mgp <- get(markers)
    }
    else{
      #if subclass_both specified run CgG markers on prefrontal and MTG on temporal
      if(cohort == "ROSMAP" | cohort =="BM10" | cohort =="BM44"){
       markers_mgp <- subclass_both_CgG_common_final
      }
      else{
        markers_mgp <- subclass_both_MTG_common_final
      }
    }
    for(cell in names(markers_mgp)){
      mgp_df[,cell] <- as.numeric(scale(mgp_df[,cell], center = TRUE, scale = TRUE))
    }
    mgp_ZScored_name <- paste0(mgp_name, "_ZScored")
    assign(mgp_ZScored_name, mgp_df)
    saveRDS(get(mgp_ZScored_name), paste0(mgp_ZScored_name, ".rds"))
    setwd('../')
  }
}

Now that we’ve Z scored all the MGPs, we can move onto our mega-analysis. Before we do though let’s use our common marker lists to determine how each of the genes that define the cell-types are behaving relative to AD diagnosis, i.e. is there a strong unidirectional association between all of the genes in a cell-type and AD diagnosis?

frenchFryPlot<-function(AD_coef, AD_pval, marker_list, cohort){
  
  
  pathology_df <- AD_pval
  
  marker_df <-  data.frame(unlist(marker_list, use.names=F),rep(names(marker_list),
                                                                lengths(marker_list)))
  colnames(marker_df) <- c("marker", "celltype")
  for (marker in names(marker_list)){
    df <- pathology_df %>% filter(marker == marker)
    if (marker == names(marker_list)[1]){
      result <- df
    }
    else{
      result <- rbind(result, df)
    }
  }
  
  final_result <- merge(result, AD_coef, by="gene")
  final_result$signedP <- -log10(final_result$pval) *(sign(final_result$coef))
  
  merge_markers<- marker_df %>% rename( gene = marker)
  final_result <- merge(merge_markers, final_result)
  final_result$cohort <- c(cohort)
  return(unique(final_result))
  
}

cohorts <- c("ROSMAP", "MAYO", "BM10", "BM22", "BM36", "BM44")
#marker_lists <- c("in_ex_CgG_common_final","subclass_CgG_common_final", "subclass_MTG_common_final","subclass_both_common_final", "lake", "darmanis")
marker_lists <- c("subclass_both_common_final")
setwd('./geneAnno')
gene_anno <- readRDS("gene_anno.rds")
setwd('../')

for (markers in marker_lists){
  for(cohort in cohorts){
    if(str_detect(markers, "subclass_both")){
      #if subclass_both specified run CgG markers on prefrontal and MTG on temporal
      if(cohort == "ROSMAP" | cohort =="BM10" | cohort =="BM44"){
       markers_for_plot <- subclass_both_CgG_common_final
      }
      else{
        markers_for_plot <- subclass_both_MTG_common_final
        }
    }
    else{
      markers_for_plot <- get(markers)
    }
    print(cohort)
    lmod_name <-  paste0("lmod_" ,cohort)
    
    filename <- paste0(lmod_name, ".rds")
    setwd('./cohortQCMods')
    lmod <- readRDS(filename)
    setwd('../')
    assign(lmod_name, lmod)
    eb <- eBayes(lmod,robust = T)
    
    if(cohort == "ROSMAP"){
      n=11
    }
    if(cohort == "MAYO"){
      n=24
    }
    if(cohort == "BM10"){
      n=18
    }
    if(cohort == "BM22"){
      n=21
    }
    else if(cohort == "BM36" || cohort == "BM44"){
      n=19
    }

    print(n)
    gene_v_AD  <-lmod$coefficients[,c(n)]
    gene_v_AD <- as.data.frame(gene_v_AD)
    gene_v_AD <- tibble:: rownames_to_column(gene_v_AD, var="Gene")
    
    final_df <- merge(gene_v_AD, gene_anno)
    final_df$new_gene<- final_df$Hgnc_Gene %>% make.names(unique = T) #this is ONE way of dealing with the duplicates, just making them into separate, unique names
    #add new_gene column w/ Hgnc_Gene values filtered to have no duplicates
    AD_coef <- final_df[,c(2,4)]
    colnames(AD_coef) <- c("coef", "gene")
    
    coef_df <- AD_coef #dataframe with coefficient column and HGNC gene name column 

    #get p values
    p_geneAD <- eb$p.value
    p_geneAD  <-p_geneAD[,c(n)]
    p_geneAD <- as.data.frame(p_geneAD)
    p_geneAD <- tibble:: rownames_to_column(p_geneAD, var="Gene")
    final_df <- merge(p_geneAD, gene_anno)
    new_gene <- final_df$Hgnc_Gene %>% make.names(unique = T) #this is ONE way of dealing with the duplicates, just making them into separate, unique names
    #add new_gene column w/ Hgnc_Gene values filtered to have no duplicates
    final_df$new_gene <- new_gene
    AD_pval <- final_df[,c(2,4)]
    colnames(AD_pval) <- c("pval", "gene")
    
    gene_pathology_association <- frenchFryPlot(AD_coef, AD_pval, cohort, marker_list =markers_for_plot)
    assign(paste0("fry_plot", cohort), gene_pathology_association)
  }
  fry_df <- rbind(fry_plotROSMAP, fry_plotMAYO, fry_plotBM10,
                  fry_plotBM22, fry_plotBM36, fry_plotBM44)
  fry_df$fdr <- p.adjust(fry_df$pval, method="fdr")
  fry_df$signedFDR <- -log10(fry_df$fdr) *(sign(fry_df$signedP))
  fry_df$sig <- ifelse(fry_df$pval < 0.05, "#94C973", "#808080")
  
  for (cell in names((markers_for_plot))){
    celltype_fry_df <- fry_df %>% filter(celltype == cell)
    celltype_heat_map <- ggplot(celltype_fry_df, aes(cohort, gene, fill= signedP))+
      theme_minimal() + geom_tile() + 
      scale_fill_gradient2(low="darkblue", high="darkgreen", guide="colorbar") +   
      ggtitle(paste0("Significance of \n" , cell, " ", markers , 
                     "\n marker genes per cohort")) 
    
    celltype_french <- ggplot(celltype_fry_df, aes(x= gene,y= signedP))+
      theme_minimal() + theme(axis.text.x = element_text(angle = 45)) +
      geom_bar(stat="identity", fill = celltype_fry_df$sig)+ 
      facet_wrap(~cohort, scales = 'free_x',nrow=6) + 
      ggtitle(paste0("Significance of \n" , cell, " ", markers , 
                     "\n marker genes per cohort"))
    
    celltype_heat_map
    celltype_french
    
    setwd('./frenchHeat')
    heat_name <- paste0("heat_map",make.names(cell))
    assign(heat_name, celltype_heat_map)
    print(get(heat_name))
    saveRDS(get(heat_name), paste0(heat_name, ".rds"))
    ggsave(paste0(heat_name, markers, "_plot", ".png"), width = 15, height = 12)
    
    french_name <- paste0("french_fry",make.names(cell))
    assign(french_name, celltype_french)
    print(get(french_name))
    saveRDS(get(french_name), paste0(french_name, ".rds"))
    ggsave(paste0(french_name, markers, "_plot", ".png"), width = 15, height = 20)
    setwd('../')
    
  }
}

## [1] "ROSMAP"
## [1] 11
## [1] "MAYO"
## [1] 24
## [1] "BM10"
## [1] 18
## [1] "BM22"
## [1] 21
## [1] "BM36"
## [1] 19
## [1] "BM44"
## [1] 19